function v = pre_pare (u, mode)

   ## usage: v = pre_norm (u, mode)
   ## 
   ## remove/add unit details from/to u.x (back and forth)
   ## FIXME: use beta distribution for most variables

   pkg load statistics optim

   global PAR clim

   cldcdf = @(p,x) betacdf(x, p(1), p(2)) ; p0 = [1 0.5]' ;

   if nargin < 2, mode = "dir" ; endif

   v = u ;

   if strcmp(mode, "inv")

      for j = 1:columns(u.x)

	 Idfd = find(!isnan(u.x(:,j))) ; v.x(:,j) = nan ;
	 y = u.x(Idfd,j) ; name = u.names{j} ; uv = u.vars{j} ;

	 switch name
	    case {"cloud" "cloudcover"}
	       a = 0 ; b = 8 ;
	       p = clim.(uv).p ;
	       y = b * betainv(normcdf(y), p(1), p(2)) ;
	       v.x(Idfd,j) = max(a, y) ;
	    case "humidity"
	       a = 0 ; b = 100 ;
	       sigma = clim.(uv).sigma ;
	       y = raylinv(normcdf(y), sigma) ;
	       v.x(Idfd,j) = min(b, max(a, b - y)) ;
	    case {"sun" "sunduration" "insolation"}
	       a = 0 ;
	       t = date2toy(u.id(Idfd,:), u.cal)*365 / (2*pi) + 10 ;
	       dl = day_length(t, u.lats(j)) ;
	       p = clim.(uv).p ; q = clim.(uv).q ;
	       y = normcdf(y) ; Ia = (y > q) ;
	       y(Ia) = (y(Ia) - q) / (1-q) ;
	       y(Ia) = betainv(y(Ia), p(1), p(2)) ;
	       y(Ia) = dl(Ia) .* y(Ia) ;
	       y(!Ia) = 0 ;
	       v.x(Idfd,j) = y ;
	    case {"wind" "windspeed"}
	       ym = clim.(uv).min ; sigma = clim.(uv).sigma ;
	       I = (y >= ym) ;
	       y = raylinv(normcdf(y(I)), sigma) ;
	       v.x(Idfd(I),j) = y ; v.x(Idfd(!I),j) = 0 ;
	    otherwise
	       v.x(:,j) = u.x(:,j) ;
	 endswitch

      endfor

   else

      for j = 1:columns(u.x)

	 Idfd = find(!isnan(u.x(:,j))) ; v.x(:,j) = nan ;
	 y = u.x(Idfd,j) ; name = u.names{j} ; uv = u.vars{j} ;

	 switch name
	    case {"cloud" "cloudcover"}
	       a = 0 ; b = 8 ;
	       [F I] = cdf(y=y/b) ;
	       [p, fy, cvg, outp] = nonlin_curvefit(cldcdf, p0, y, F) ; clim.(uv).p = p ;
	       y = betacdf(y, p(1), p(2)) ; Ia = (0 < y) ; Ib = (y < 1) ;
	       y = norminv(y(Ia & Ib)) ; ya = min(y) ; yb = max(y) ;
	       v.x(Idfd(Ia & Ib),j) = y ; v.x(Idfd(!Ia),j) = ya ; v.x(Idfd(!Ib),j) = yb ;
	    case "humidity"
	       a = 0 ; b = 100 ;
	       I = (y < b) ;
	       sigma = mean(y=b-y(I)) * sqrt(2/pi) ; clim.(uv).sigma = sigma ;
	       y = norminv(raylcdf(y, sigma)) ;
	       v.x(Idfd(I),j) = y ; v.x(Idfd(!I),j) = min(y) ;
	    case {"sun" "sunduration" "insolation"}
	       a = 0 ;
	       t = date2toy(u.id(Idfd,:), u.cal)*365 / (2*pi) + 10 ;
	       dl = day_length(t, u.lats(j)) ;
	       I = (y > a) ; q = sum(y<=a) ; q = clim.(uv).q = q / (q + sum(I)) ;
	       [F II] = cdf(y=y(I)./dl(I)) ;
	       [p, fy, cvg, outp] = nonlin_curvefit(cldcdf, p0, y, F) ; clim.(uv).p = p ;
	       y = q + (1 - q) * cldcdf(p, y) ; Ib = (y < 1) ;
	       y(Ib) = norminv(y(Ib)) ; y(!Ib) = max(y(Ib)) ;
	       v.x(Idfd(!I),j) = trunc_rnd("norm", sum(!I), -Inf, norminv(q)) ;
	       v.x(Idfd(I),j) = y ;
	    case {"wind" "windspeed"}
	       a = 0 ; b = inf ;
	       sigma = mean(y) * sqrt(2/pi) ; clim.(uv).sigma = sigma ;
	       I = (y > a) ;
	       y = norminv(raylcdf(y(I), sigma)) ; clim.(uv).min = ym = min(y) ;
	       v.x(Idfd(I),j) = y ; v.x(Idfd(!I),j) = ym ;
	    otherwise
	       v.x(:,j) = u.x(:,j) ;
	 endswitch

      endfor

   endif

endfunction
